Modulational instability in Bose-Einstein condensates in optical lattices 
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A self consistent theory of a cylindrically shaped Bose-Einstein condensate (BEC) periodically mod- 
ulated by a laser beam is presented. We show, both analytically and numerically, that modulational 
instability/stability is the mechanism by which wavefunctions of soliton type can be generated in 
cylindrically shaped BEC subject to a one-dimensional optical lattice. The theory explains why 
bright solitons can exist in BEC with positive scattering length and why condensate with negative 
scattering length can be stable and give rise to dark solitary pulses. PACS numbers: 42.65.-k, 42.50. 
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There is an increasing interest in the study of Bose- 
Einstein condensates (BEC) in presence of periodic po- 
tentials, such as the one induced by detuned standing 
waves of light (optical lattices) Q. Switching on an opti- 
cal lattice in a continuous BEC induces fragmentation of 
the original wavefunction into local wavefunctions cen- 
tered around the minima of the potential, this leading 
to a crystal-like structure of mutually interacting BECs. 
In analogy with the usual theory of crystals, one can 
think to control the dynamics of this new state of mat- 
ter by properly choosing the parameters of the lattice. 
This gives, for example, the possibility to observe macro- 
scopic quantum interference phenomena with emission of 
coherent pulses of atoms (Bloch oscillations), as recently 
reported in Ref. || for vertical BEC arrays in the grav- 
itational field. Understanding the properties of BEC in 
optical lattices is therefore of fundamental importance 
for developing novel applications of quantum mechanics 
such as atom lasers and atom interferometers. For small 
overlapping between local wavefunctions, a tight-binding 
model can be developed. This was done, for the one di- 
mensional (ID) case, in Refs. || where it was shown that 
the mean field equation for the condensate wavefunction 
reduces to the so called discrete nonlinear Schrodinger 
equation Q|. The tight-binding approximation, however, 
putting restrictions on the shape of the wave function 
(i.e. on the number of atoms in the condensate), as well 
as on the potential profile, is applicable only to partic- 
ular experimental settings. From this point of view it 
is desirable to develop a theory of BEC in optical lat- 
tices which does not rely on this approximation. Studies 
in this direction were made in terms of a ID nonlinear 
Schrodinger equation (NLS) with trigonometric Q or el- 
liptical potentials || . Bright and dark solitons in BEC in 
optical lattices, analogue to the gap-soliton of photonic 
crystals , were also shown to exist H|] . 

The aim of this Letter is to investigate, both analyti- 
cally and numerically, modulational instability phenom- 
ena of extended states at the border of the Brillouin zone. 
To this end we construct approximate ground state so- 
lutions of the original 3D by means of a multiple scale 



expansion, starting from the exact eigenfunctions of the 
underlying linear Schrodinger equation with potentials 
which are parabolic in the transverse direction and pe- 
riodic in the longitudinal one (periodic cylindrical trap). 
We show that at the lowest orders in the expansion the 
condensate evolves according to an effective ID NLS with 
the dispersive term depending on the effective mass of the 
Bloch states of the underlying linear problem. Extended 
states close to the borders of the Brillouin zone, are then 
shown to be unstable (stable) against small spatial mod- 
ulations (modulational instability) depending on the sign 
of the dispersion in the effective ID NLS. The stability 
properties of these states is shown to be the basic mech- 
anism by which bright (resp. dark) solitons are created 
in BEC with positive (resp. negative) scattering lengths. 
Numerical simulations of the longitudinal BEC dynamics 
confirm the predictions of our theory. The possibility to 
observe the modulational instability phenomena in real 
BEC is discussed at the end of the paper. 

As is well known Jicj ], the condensate wavefunction is 
described by the Gross-Pitaevskii equation (CPE) 
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with go = 4iTh 2 a s /m, m is the atomic mass, and a s is 
the s-wave scattering length of atoms which can be ei- 
ther positive or negative. We consider a trap potential 
of the form V(r) = im^ 2 r^_ + Vb cos(rez), which model a 
cylindrically shaped BEC periodically modulated along 
the z-axis (the results, however, will not depend on 
the form of periodic potential used, and can be eas- 
ily generalized to arbitrary z-periodic potentials). Here 
r = (r±,z), Vb is the potential deepness, v the trap fre- 
quency in the transverse direction, and 2n/n the period 
of the modulation. We assume periodic boundary condi- 
tions ^(r±, z,t) = i f?(r±, z + L,t), with L denoting the 
length of the cylinder. The change of variables 1 1— > 2t/v, 
r i ► a r, * h-> (iV/aj]) 1 / 2 ^, with a = (h/ '(mi/)) 1 / 2 , al- 
lows us to rewrite Eq.(Q) in the dimensionless form 
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where x = 8nNa s /ao, and C = C± + C z with 



-d 2 /dz 2 + 2Acos(fcz), (3) 



(here A_l denotes the two-dimensional Laplacian, k — 
ao/n and A = Vqu/%). In these units the wave function 
results normalized to one, i.e. 



dr± 



dz[i,[ 2 = 1, 



(4) 



with L = L/ao denoting the normalized length of the 
cylinder. In the following we shall restrict to the small 
amplitude limit (xlV'l 2 ^ 1) an d construct a solution of 
Eq. (H) perturbatively, starting from the solution of the 
linear problem. These last can be written as products of 
eigenfunctions of the operators in Eq. (0) 



£>z<f>nq(z) = £nq<j>nq{z), £j_£ nm (rjj 



For the considered potential, <pa q (z) are solutions of the 
Mathieu equation, while £nm(r_i_) are eigenfunctions of 
the two-dimensional harmonic oscillator (n and m denote 
the principal and the angular quantum numbers of the 
harmonic oscillator, while h and q denote the band index 
and the wavevector inside the first Brillouin zone of the 
ID lattice, respectively). We look for solutions of Eq. (||) 
of the form 



1> = 



(aipi + a 2 4>2 H ), 



(5) 



with a a small parameter whose physical meaning will 
be clarified later (the prefactor is unimportant and in- 
troduced just for convenience). Since we are interested 
in the ground state we take as the leading order term in 
Eq. (||) a small modulation of the linear ground state 
wavefunction {no = 0, mo = 0, no = 1) of the form 
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Lo{q). The modu- 



with u). 

lating amplitude A(z, t) is considered to be a function 
of a set of independent spatial and temporal variables 
of the form z = (z\, z%,..., z n , . . .) with z n = a n z, and 
t = (ti,tz, . . . ,t n , . . .) with t n = a n t, respectively. To 
simplify the notation we introduce the shortcut symbols 

£0 = £n ,m , £{q) = £n ,q, 4>q( Z ) = 4>n ,q(z), and in 

the modulation amplitude A, we show only the depen- 
dence on the most "rapid" variables. The time and co- 
ordinate derivatives in Eq. (^J) are then expanded as 
d/dt = Y. a =^ a d/dt a and 3/dz = E a ^"9/&«. 
Substituting the above expansions in Eq. (0) and col- 
lecting all the terms of the same order in a, we obtain at 
the first order: idipi/dto — Ci/j\ — 0, which is evidently 
satisfied by ip\ given by (||). At the second order in <r, 
the following equation is obtained 



,dip2 „ .dipi d 2 tpi 



(7) 



dto dti ~dz dzi' 

whose solution can be searched in the form 

n,m (h,q')^(h ,q) 

Substituting (||) in (0) and projecting along the eigen- 
functions of operators (||) with h ^ ho, we find that 



-iuj(q)t 



(8) 
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with T fvho (q) = -2 f Q <j>nq{z)-^4>n oq { z ) dz. The solvabil- 
ity condition of Eq. Q reads as |^ 
which we see that A = A(^; z%, £2); with £ = Z\ 
Note that v = v(q) = iTn n (q) can be interpreted as 
the group velocity of the wave packet in the z-direction. 
Finally, at the third order in a, we get 



u|4. = 0, from 
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, d 2 j> 2 
dz dzi 



(10) 



Requiring orthogonality (to avoid secular terms) between 
the right hand side of this equation and the kernel of 
the operator id/dto — C, and taking into account the 
expressions of ipi and tp2 derived above, we find that 
Eq.(|lC)|) reduces to the following NLS equation 
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DSi+ x [A{ 2 A = 0, (11) 
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where D = D(q) = 1 + V- 

fective group velocity dispersion induced by the periodic 
potential, and 
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is the effective nonlinearity (here we integrated on ra- 
dial variables and used the ground state wavefunction of 
the 2D harmonic oscillator). The above expressions of 
v and D, in terms of eigenfunctions of the linear opera- 
tor £, can be simplified by expressing them in terms of 
the energy spectrum of the non interacting linear system. 
This can be done in the same manner as in the theory of 
optical gap solitons Qj. To this end, we take two close 
Bloch solutions of the ID linear problem, of the form 
4>q{z) = cxp(iqz)un_ q (z) , which differ only by a small Sq, 
so that Un, q +8q{z) can be considered as a perturbation of 
Un. q (z) generated by the operator —2iSq (-4- + iq) + {5q) 2 . 
This perturbation produces a shift A = £n, q +s q — £n. q 
in energy, which can be expanded in a Taylor series in 
Sq. On the other hand, A can also be computed from 
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perturbation theory. A comparison of the correspond- 
ing expressions leads to v = duj(q)/dq and D = | d £^ , 
i.e. v and D are, respectively, the slope (velocity) and 
the curvature (inverse effective mass) of the energy band 
(Bloch states) of the underlying linear problem. 

From the physical point of view the above results have 
a number of consequences. First, the group velocity in- 
duced by the periodicity at the boundaries of zone dom- 
inate the dispersion inherent to NLS. For example, if we 
take k = 2.0 and A = 0.5, we have that the edges of the 
first gap £( 2 )] are at £W ~ .47, and £^ « 1.47. 
The effective dispersion at these points is oj" ~ —6.13, 
and u>2 w 10.14, respectively (here oj'J — d 2 uj/dq 2 \ q=qj ). 
Thus even in the case the group velocity dispersion does 
not change sign it becomes much larger than the NLS dis- 
persion. Second, for fixed nonlinearity and in presence 
of the periodic potential, the dynamics will crucially de- 
pend on the sign of D. This sign can be controlled by 
changing the wave number of the initial state, as well as, 
the potential parameters. Instability phenomena of ex- 
tended (Bloch) states close to the edges of the Brillouin 
zone can then appear. To understand this, let us assume 
positive scattering length (x > in (f[l])) and consider 
the Bloch state at S^ 1 ', for which £)W < 0. In the pres- 
ence of a repulsive inter-atomic interaction (% > 0), the 
energy of this state will be shifted upward in the gap 
where it cannot exists. One can expect then the state to 
become unstable against small spatial modulations (mod- 
ulational instability) so that new excitations must arise. 
Equation (pd] ) predicts that out of the instability bright 
solitons should appear (recall that for \ > 0, and D < 
(resp.D > 0), Eq. ( |ll| ) has stable bright (resp. dark) soli- 
ton solutions) . On the contrary, if we take as initial state 
the Bloch state at the bottom of the second band, £^ 2 \ 
where D = > 0, one expects modulational stability 
instead (in this case the nonlinearity is pulling the energy 
of the state further up in the second band where it can 
still exist). This extended stable state can be then used 
as background to construct the dark soliton solution ex- 
pected in this case from Eq.(pd|) (see below). Obviously, 
for negative scattering lengths the opposite situation will 
occur i.e. modulational instability will appear at the top 
of the gap (leading to bright solitons) and stability at 
the bottom (leading to dark solitons). From this it is 
clear that the stability properties of the Bloch states at 
the edge of the Brillouin zone, plays a crucial role for the 
existence of bright and dark solitons in BEC in optical 
lattices both for positive and negative scattering lengths. 

These predictions can be easily checked by direct nu- 
merical integration. To this regard we remark that insta- 
bilities along z-direction mainly depend on the spectrum 
of the operator C z (the transverse distribution of the con- 
densate affects only the absolute value of the coefficient 
X), so that we can perform numerical simulations in the 
framework of a ID NLS equation obtained from Eq. (||) 



with L ~ C z . Moreover, we note that the Bloch state 
(Mathieu function) at the top of the first band (bottom 
of the gap), is an odd function of z which can be ap- 
proximated by sin(z), while the one at the bottom of the 
second band (top of the gap) is an even function of z 
very close to cos(z). In the following we shall use these 
approximate states as initial conditions for investigating 
modulational stability since they are, in real experiments, 
more easy to generate. 

FIG. 1. Modulational instability in Eq. (§]) with £ « C z 
for parameter values x — L0, k = 2.0 and A = 0.5. The 
initial condition is an approximated eigenfunction, taken as 
a sine function, of the first band of the linear system at the 
edge £^ w 0.47 of the Brillouin zone. 

In Fig. [l] a numerical simulation of the ID problem 
with initial condition close to the state at the bottom 
of the gap, is depicted. We see that, as expected from 
our analysis, modulational instability develops and, in 
spite of the fact that we have positive scattering (x = 1), 
bright solitons are created in agreement with our analy- 
sis (the number of solitons coming out from the instabil- 
ity can be estimated as Lk max /(2Tr), where k max is the 
wavenumber of the most unstable linear mode (Tip). We 
remark that although the theory is valid for small am- 
plitude excitations, the numerical simulations show that 
the obtained results extend also above this limit (note 
that in Fig. |l| x = !)• An intuitive explanation for this 
is that small amplitude solitons once formed can only 
become more and more localized as the nonlinearity is 
increased. The modulational instability at higher nonlin- 
earity should therefore produce solitons which are more 
localized and of large amplitude. This is precisely what 
observed in Fig 1. In contrast to this, we find that an 
initial condition corresponding to a Bloch state close to 
the top of the gap, remains modulationally stable also 
in presence of nonlinearity. This is reported in Fig. ^ 
for an initial profile of cosine type. It is interesting to 
note that one can use this state to construct the sta- 
ble dark soliton predicted by Eq. ([ll]). To this end we 
take as initial condition a modulated Bloch state of the 
form tanh(Az) cos(z), where the cosine function, taken as 
background, approximates the Mathieu eigenfunction at 
the edge E^ 2 \ while the tanh modulation is used to make 
the profile close to the expected dark state. 

FIG. 2. Same as in Fig.l but for the eigenfunction at the 
top of the gap £^ ~ 1.47 approximated with a cosine func- 
tion. 

FIG. 3. Same as in Fig.l but for a dark soliton initial con- 
dition 

In Fig. |3| the corresponding numerical simulation is re- 
ported, from which we see that a dark soliton is indeed 
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generated, in perfect agreement with our analysis (the 
energy of this state is in the gap close to the bottom of 
the second band). It is interesting to note that for nega- 
tive scattering lengths this leads to the existence of dark 
soliton in BEC in optical lattices (in this case one must 
use the stable Bloch state at the bottom of the gap as 
background for the dark solution). 

In order to check the self-consistency of the theory we 
shall estimate the size of the parameter a used for the ex- 
pansion, and the magnitude of the effective nonlinearity 
in Eq. (|ll|) . To this end we start with the dark soliton or 
periodic solution and notice that the eigenfunctions <j> q {z) 
are normalized to one, so that L\4> q \ 2 ~ 1 and hence, from 
Eq. ( |l2| ) we have that x ~ 1. Similarly, from the nor- 
malization of the wavefunction (Q) and from the expan- 
sion (||), we have that nL 2 a 2 /\x\ ~ 1 from which, after 
restoring physical units, we get a 2 = 8Na s ao/L 2 . If we 
consider the case of a condensate with N « 10 4 atoms 
of 87 Rb (a a ~ 5.5nm) with a radial size ao ~ 17 /zm, and 
length L ~ 300 fim we have that a 2 ~ 0.08, this being 
reasonably small to justify our expansion (smaller values 
can be achieved by considering longer, thinner, cylin- 
ders and smaller values of N). As to the initial state, 
we remark that it could be generated from an uniform 
cylinder by modulating it along the z-axis with a sine 
or cosine wave of light with twice the wavelength of 
the lattice. Another possibility is to use an initial uni- 
form condensate and accelerate the lattice until the state 
reaches the edge of the band (it is enough to be close 
to the edge for the instability to develop). Finally we 
mention, that although the expansion has been provided 
for a cylindrically-shaped BEC, a number of effects dis- 
cussed is relevant to a cigar-shaped BEC (i.e. including 
parabolic confining potential in the direction of periodic- 
ity). This is the case when the effect (instability, bright 
or static dark soliton) observed has a scale much less than 
the length of the condensate. We hope that the phenom- 
ena of modulational instability discussed in this Letter 
will be soon observed in real BEC experiments. 
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